%load_ext autoreload
%autoreload 2
import math
import additive2.utility as u
import additive2.clusters as cl
import matplotlib.pyplot as plt
import joblib  # type: ignore
import pandas as pd
import numpy as np
import seaborn as sns
u.custom_matplotlib_style()
from IPython.display import display  # type: ignore
import additive2.notebooks.pore_preformance as nb

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

Pore performance

Now let's see how the detected pores differ in terms of pixels. That is if a pixel is detected as pore pixel in the high-res, is it also detected in low res?

# break point
pores_mid_res = joblib.load("../data/Additive/pores-mid-res-v01.np")
pores_high_res = joblib.load("../data/Additive/pores-high-res-v01.np")
pores_pred = joblib.load("../data/Additive/pores-pred-v01.np")
mask = joblib.load("../data/Additive/mask-v01.np")

# let's try to put the three in a single plot
fig = plt.figure(figsize=plt.figaspect(0.3))
axes = []
for i in range(1, 4):
    axes.append(fig.add_subplot(1, 3, i, projection='3d'))

pore_imgs = [pores_mid_res, pores_pred, pores_high_res]
titles = ['Mid-resolution', 'Improved resolution', 'High-resolution']


for ax, pore_img, title in zip(axes, pore_imgs, titles):
    u.draw_pores(pore_img, mask, ax=ax)
    ax.set_title(title)
    ax.axis('off')

Pore count

let's find a way to compute pairwise pixel distance

clusters_high_res, indices_high_res = cl.get_clusters(pores_high_res)
clusters_mid_res, indices_mid_res = cl.get_clusters(pores_mid_res)
clusters_pred, indices_predicted = cl.get_clusters(pores_pred)

indices_high_res.shape

(1711, 3)

clusters_mid, pts_mid, clusters_mid_ = cl.get_pores(pores_mid_res)

clusters_mid_res.max(), clusters_pred.max(), clusters_high_res.max()

(1390, 1865, 3013)

clusters_pred, pts_pred, clusters_pred_ = cl.get_pores(pores_pred)

clusters_high, pts_high, clusters_high_ = cl.get_pores(pores_high_res)

assert len(clusters_mid) == len(pts_mid)
assert len(clusters_high) == len(pts_high)
assert len(clusters_pred) == len(pts_pred)

Matching Pores

pred_high_map, high_pred_map = cl.find_matching_pores(pts_pred, pts_high)
mid_high_map, high_mid_map = cl.find_matching_pores(pts_mid, pts_high)

# the following should match pixel-wise comparisons above
(pred_high_map>-1).sum() / len(pred_high_map), (high_pred_map>-1).sum() / len(high_pred_map), (mid_high_map>-1).sum() / len(mid_high_map), (high_mid_map>-1).sum() / len(high_mid_map)

(0.7715996578272027, 0.527177089421391, 0.5581151832460733, 0.3115137346580947)

# the following should match pixel-wise comparisons above
def get_precision_recall(low, high):
    tp = np.sum(low>-1)
    fp = np.sum(low==-1)
    fn = np.sum(high==-1)
    return tp/(tp+fp), tp/(tp+fn)


precision_recall_df = pd.DataFrame(
    [get_precision_recall(pred_high_map, high_pred_map), get_precision_recall(mid_high_map, high_mid_map)],
    columns = ['Precision', 'Recall']
)
precision_recall_df.index = ['Improved Resolution', 'Medium Resolution']
precision_recall_df['Accuracy'] = 1
precision_recall_df.round(2)

	Precision	Recall	Accuracy
Improved Resolution	0.77	0.53	1
Medium Resolution	0.56	0.31	1

Which it does.

Now let's compare the three images based on the numbers of matching pores

cl_pred_high_map = cl.find_matching_clusters(pred_high_map, clusters_pred, clusters_high)
cl_mid_high_map = cl.find_matching_clusters(mid_high_map, clusters_mid, clusters_high)

# total number of matched pores between mid-res and high res
len(cl_mid_high_map)

14

# total number of matched pores between predicted and high res
len(cl_pred_high_map)
# let's say it is 25

18

# total number of pores detected by high-res is:
np.max(clusters_high)

38

let's check how individual pore sizes compare to each other

# comparing mid vs high pores sizes
res = u.get_pore_size_diff(cl_mid_high_map, clusters_mid_, clusters_high_, 
                           name1='mid', name2='high')
display(res)
res['rel. error'].describe()

	mid	high	rel. error
0	47	54	0.15
1	36	72	1.00
2	39	51	0.31
3	47	61	0.30
4	56	57	0.02
5	37	37	0.00
6	35	61	0.74
7	45	78	0.73
8	44	68	0.55
9	63	79	0.25
10	25	45	0.80
11	33	76	1.30
12	23	46	1.00
13	57	60	0.05

count    14.000000
mean      0.514286
std       0.419170
min       0.000000
25%       0.175000
50%       0.430000
75%       0.785000
max       1.300000
Name: rel. error, dtype: float64

res = u.get_pore_size_diff(cl_pred_high_map, clusters_pred_, 
                           clusters_high_, name1='pred', name2='high')
display(res)
res['rel. error'].describe()

	pred	high	rel. error
0	52	55	0.06
1	56	54	0.04
2	66	72	0.09
3	53	51	0.04
4	63	61	0.03
5	37	43	0.16
6	56	57	0.02
7	39	37	0.05
8	50	61	0.22
9	73	78	0.07
10	67	68	0.01
11	42	61	0.45
12	28	81	1.89
13	75	79	0.05
14	36	45	0.25
15	63	76	0.21
16	38	46	0.21
17	55	60	0.09

count    18.000000
mean      0.218889
std       0.431903
min       0.010000
25%       0.042500
50%       0.080000
75%       0.210000
max       1.890000
Name: rel. error, dtype: float64

Getting pore features

# pts_mid
pts_ = pts_mid[np.array(clusters_mid) == 0]

pf = u.PoreFeatures(pts_)

mid_features, high_features_mid = u.get_feature_pairs(cl_mid_high_map, pts_mid, 
                                                      clusters_mid, 
                                                    pts_high, clusters_high)
pred_features, high_features_pred = u.get_feature_pairs(cl_pred_high_map, pts_pred, clusters_pred, 
                                                 pts_high, clusters_high)

u.ttest_1samp([1, 2, 3, 4, 5, 6], 0)

Ttest_1sampResult(statistic=4.58257569495584, pvalue=0.00593354451759226)

stats_ = pd.DataFrame(
    (high_features_pred.values - pred_features.values),
    columns=pred_features.columns
).assign(tmp=1).groupby('tmp').agg(['mean', 'count', 'std', 
                                    u.ttest_pvalue]).T.reset_index()
stats = pd.pivot_table(stats_, index='level_0', columns='level_1', 
                       values=1).drop('cluster')
ci95_hi = []
ci95_lo = []
for i in stats.index:
    m, c, s, *_ = stats[['mean', 'count', 'std']].loc[i]
    ci95_hi.append(m + 1.96*s/math.sqrt(c))
    ci95_lo.append(m - 1.96*s/math.sqrt(c))

stats['ci95_lo'] = ci95_lo
stats['ci95_hi'] = ci95_hi
display(stats.round(2).sort_index())

level_1	count	mean	std	ttest_pvalue	ci95_lo	ci95_hi
level_0
Convex Hull Area	18.0	8.58	18.07	0.06	0.23	16.92
Convex Hull Volume	18.0	6.38	17.03	0.13	-1.49	14.25
Ellipsoid Area	18.0	8.47	25.25	0.17	-3.20	20.14
Ellipsoid Elongation	18.0	-0.05	0.84	0.82	-0.44	0.34
Ellipsoid Volume	18.0	10.46	35.26	0.23	-5.83	26.75
Pore Elongation	18.0	0.03	0.94	0.89	-0.40	0.47
Pore Sphericity	18.0	0.00	0.04	0.59	-0.01	0.02
Pore Volume	18.0	12.72	23.33	0.03	1.94	23.50

stats_ = pd.DataFrame(
    (high_features_mid.values - mid_features.values),
    columns=mid_features.columns
).assign(tmp=1).groupby('tmp').agg(['mean', 'count', 'std', 
                                    u.ttest_pvalue]).T.reset_index()
stats = pd.pivot_table(stats_, index='level_0', columns='level_1', values=1).drop('cluster')
ci95_hi = []
ci95_lo = []
for i in stats.index:
    m, c, s, *_ = stats[['mean', 'count', 'std']].loc[i]
    ci95_hi.append(m + 1.96*s/math.sqrt(c))
    ci95_lo.append(m - 1.96*s/math.sqrt(c))
    

stats['ci95_lo'] = ci95_lo
stats['ci95_hi'] = ci95_hi
display(stats.round(2).sort_index())

level_1	count	mean	std	ttest_pvalue	ci95_lo	ci95_hi
level_0
Convex Hull Area	14.0	32.91	39.55	0.01	12.20	53.63
Convex Hull Volume	14.0	14.86	17.65	0.01	5.61	24.10
Ellipsoid Area	14.0	37.26	60.40	0.04	5.62	68.90
Ellipsoid Elongation	14.0	1.02	1.67	0.04	0.15	1.90
Ellipsoid Volume	14.0	21.08	34.09	0.04	3.22	38.93
Pore Elongation	14.0	1.36	1.70	0.01	0.47	2.25
Pore Sphericity	14.0	-0.02	0.03	0.01	-0.04	-0.01
Pore Volume	14.0	38.14	32.31	0.00	21.22	55.07

# todo write about the model

plotting individual matched pores

import additive2.plotEllipsoid as et
plt.figure(figsize=(15, 15))
fig = plt.figure(figsize=(24, 8))
axes = []
for i in range(1, 4):
    axes.append(fig.add_subplot(1, 3, i, projection='3d'))
    
ax = axes[0]
pf = u.PoreFeatures(pts_mid[np.array(clusters_mid)==27])
ellipsoid_tool = et.EllipsoidTool()
u.draw_voxels(pf.pts_in, ax=ax, alpha=.1)
ellipsoid_tool.plotEllipsoid(pf.ellipsoid, plotAxes=True, ax=ax)
   
ax = axes[2]
pf = u.PoreFeatures(pts_high[np.array(clusters_high)==31])
u.draw_voxels(pf.pts_in, ax=ax, alpha=.1)
ellipsoid_tool = et.EllipsoidTool()
ellipsoid_tool.plotEllipsoid(pf.ellipsoid, plotAxes=True, ax=ax)

ax = axes[1]
pf = u.PoreFeatures(pts_pred[np.array(clusters_pred)==28])
ellipsoid_tool = et.EllipsoidTool()
u.draw_voxels(pf.pts_in, ax=ax, alpha=.1)
ellipsoid_tool.plotEllipsoid(pf.ellipsoid, plotAxes=True, ax=ax)


for ax in axes:
    ax.set_xlim([392, 406])
    ax.set_ylim([418, 432])
    ax.set_zlim([222, 233])

titles = ['Mid-resolution', 'Improved resolution', 'High-resolution']
for ax, title in zip(axes, titles):
    ax.set_title(title)

<Figure size 1080x1080 with 0 Axes>

predHighPoreMap_ = pd.DataFrame([list(x) for x in cl_pred_high_map], columns=['pred', 'high'])
midHighPoreMap = pd.DataFrame([list(x) for x in cl_mid_high_map], columns=['mid', 'high'])
predHighPoreMap = pd.concat([
    predHighPoreMap_,
    pd.DataFrame(set(range(len(clusters_high_))) - set(predHighPoreMap_['high']), columns=['high']),
    pd.DataFrame(set(range(len(clusters_pred_))) - set(predHighPoreMap_['pred']), columns=['pred']),
])

common_pores = predHighPoreMap.merge(midHighPoreMap, how='left')
common_pores = pd.concat([
    common_pores,
    pd.DataFrame(set(range(len(clusters_mid_))) - set(common_pores['mid']), columns=['mid']),
])
common_pores

	pred	high	mid
0	0.0	0.0	NaN
1	2.0	2.0	1.0
2	3.0	9.0	7.0
3	5.0	10.0	8.0
4	6.0	12.0	10.0
...	...	...	...
16	NaN	NaN	26.0
17	NaN	NaN	28.0
18	NaN	NaN	31.0
19	NaN	NaN	32.0
20	NaN	NaN	33.0

74 rows × 3 columns

import additive2.plotEllipsoid as et
plt.figure(figsize=(15, 15))
fig = plt.figure(figsize=(24, 8))
axes = []
for i in range(1, 4):
    axes.append(fig.add_subplot(1, 3, i, projection='3d'))

pred_cluster_num, high_cluster_num, mid_cluster_num = common_pores.loc[14]
    
ax = axes[0]
pf = u.PoreFeatures(pts_mid[clusters_mid_[mid_cluster_num]])
ellipsoid_tool = et.EllipsoidTool()
u.draw_voxels(pf.pts_in, ax=ax, alpha=.1)
ellipsoid_tool.plotEllipsoid(pf.ellipsoid, plotAxes=True, ax=ax)
   
ax = axes[2]
pf = u.PoreFeatures(pts_high[clusters_high_[high_cluster_num]])
u.draw_voxels(pf.pts_in, ax=ax, alpha=.1)
ellipsoid_tool = et.EllipsoidTool()
ellipsoid_tool.plotEllipsoid(pf.ellipsoid, plotAxes=True, ax=ax)

ax = axes[1]
pf = u.PoreFeatures(pts_pred[clusters_pred_[pred_cluster_num]])
ellipsoid_tool = et.EllipsoidTool()
u.draw_voxels(pf.pts_in, ax=ax, alpha=.1)
ellipsoid_tool.plotEllipsoid(pf.ellipsoid, plotAxes=True, ax=ax)


for ax in axes:
    ax.set_xlim([392, 406])
    ax.set_ylim([418, 432])
    ax.set_zlim([222, 233])

titles = ['Mid-resolution', 'Improved resolution', 'High-resolution']
for ax, title in zip(axes, titles):
    ax.set_title(title)

<Figure size 1080x1080 with 0 Axes>

import additive2.plotEllipsoid as et
dTmp = common_pores.dropna()
for j in dTmp.index:
    fig = plt.figure(figsize=(24, 8))
    axes = []
    for i in range(1, 4):
        axes.append(fig.add_subplot(1, 3, i, projection='3d'))

    pred_cluster_num, high_cluster_num, mid_cluster_num = dTmp.loc[j]

    ax = axes[0]
    pf_mid = u.PoreFeatures(pts_mid[np.array(clusters_mid)==mid_cluster_num])
    ellipsoid_tool = et.EllipsoidTool()
    u.draw_voxels(pf_mid.pts_in, ax=ax, alpha=.1)
    ellipsoid_tool.plotEllipsoid(pf_mid.ellipsoid, plotAxes=True, ax=ax)

    ax = axes[2]
    pf_high = u.PoreFeatures(pts_high[np.array(clusters_high)==high_cluster_num])
    u.draw_voxels(pf_high.pts_in, ax=ax, alpha=.1)
    ellipsoid_tool = et.EllipsoidTool()
    ellipsoid_tool.plotEllipsoid(pf_high.ellipsoid, plotAxes=True, ax=ax)

    ax = axes[1]
    pf_pred = u.PoreFeatures(pts_pred[np.array(clusters_pred)==pred_cluster_num])
    ellipsoid_tool = et.EllipsoidTool()
    u.draw_voxels(pf_pred.pts_in, ax=ax, alpha=.2)
    ellipsoid_tool.plotEllipsoid(pf_pred.ellipsoid, plotAxes=True, ax=ax)

    pf_pts = np.concatenate([pf_pred.pts, pf_high.pts, pf_mid.pts])
    mn = pf_pts.min(axis=0) - 2
    mx = pf_pts.max(axis=0) + 2

    for ax in axes:
        ax.set_xlim([mn[0], mx[0]])
        ax.set_ylim([mn[1], mx[1]])
        ax.set_zlim([mn[2], mx[2]])
        ax.set_xlabel('x')
        ax.set_ylabel('y')


    titles = ['Mid-resolution', 'Improved resolution', 'High-resolution']
    for ax, title in zip(axes, titles):
        ax.set_title(title)
    plt.show()

fig, axes = plt.subplots(1, 3, figsize=(15, 5))
axes = axes.reshape(-1)
axes[0].hist([len(x) for x in clusters_mid_], bins=5)
axes[0].set_title('Mid-resolution')
axes[1].hist([len(x) for x in clusters_high_], bins=5)
axes[1].set_title('high-resolution')
axes[2].hist([len(x) for x in clusters_pred_], bins=5);
axes[2].set_title('super-resolution')
for ax in axes:
    ax.set_xlim(10, 110)

midPoreFeatures  = pd.DataFrame([u.PoreFeatures(pts_mid[x]).all_features for x in clusters_mid_])\
    .reset_index().rename(columns={'index': 'pore'})
# midPoreFeatures.columns = ['Mid ' + x for x in midPoreFeatures.columns]

highPoreFeatures = pd.DataFrame([u.PoreFeatures(pts_high[x]).all_features for x in clusters_high_])\
    .reset_index().rename(columns={'index': 'pore'})
# highPoreFeatures.columns = ['High ' + x for x in highPoreFeatures.columns]

predPoreFeatures = pd.DataFrame([u.PoreFeatures(pts_pred[x]).all_features for x in clusters_pred_])\
    .reset_index().rename(columns={'index': 'pore'})
# predPoreFeatures.columns = ['Super ' + x for x in predPoreFeatures.columns]

len(clusters_mid_),  len(clusters_high_), len(clusters_pred_)

(35, 39, 32)

commonPoresUnstacked = common_pores.unstack().rename('pore').reset_index()\
    .rename(columns={'level_0': 'type', 'level_1': 'main_pore'})
commonPoresUnstacked.head(4)

	type	main_pore	pore
0	pred	0	0.0
1	pred	1	2.0
2	pred	2	3.0
3	pred	3	5.0

commonPoresUnstacked.groupby('type').max().sum()

main_pore    156.0
pore         103.0
dtype: float64

allPoreFeatures = pd.concat([
    midPoreFeatures.assign(type='mid'),
    highPoreFeatures.assign(type='high'),
    predPoreFeatures.assign(type='pred')
])

allPoreFeatures.groupby('type')['pore'].max()

type
high    38
mid     34
pred    31
Name: pore, dtype: int64

allSharedPoresFeatures = allPoreFeatures.merge(commonPoresUnstacked, on=['type', 'pore'])
allSharedPoresFeatures

	pore	Convex Hull Area	Pore Volume	Convex Hull Volume	Ellipsoid Volume	Ellipsoid Area	Ellipsoid Elongation	Pore Sphericity	Pore Elongation	type	main_pore
0	0	82.498975	62	30.5	71.452453	136.199815	4.375653	0.572228	5.336956	mid	0
1	1	135.113122	128	54.0	110.460570	204.864846	6.230679	0.511347	7.479158	mid	1
2	2	81.467619	70	30.0	71.158961	135.749488	4.368205	0.573121	5.312026	mid	1
3	3	77.659737	66	28.0	64.064661	124.088692	4.116301	0.574196	5.282254	mid	2
4	4	97.052389	74	36.5	99.735116	185.903807	5.711054	0.548285	6.067092	mid	3
...	...	...	...	...	...	...	...	...	...	...	...
101	27	99.680190	82	37.0	91.263904	172.922892	5.489174	0.538695	6.396917	pred	52
102	28	107.912687	104	41.5	89.201645	168.122077	5.280513	0.537168	6.451623	pred	14
103	29	208.246167	176	86.0	173.643681	317.571353	9.391903	0.452452	10.796470	pred	15
104	30	119.685728	110	47.0	97.843557	182.623290	5.625377	0.526228	6.862447	pred	16
105	31	145.714777	144	58.5	126.855739	233.815149	7.023366	0.500132	7.993679	pred	17

106 rows × 11 columns

plt.subplots(figsize=(16, 8))
sns.barplot(data=allSharedPoresFeatures, 
            x='main_pore', y='Pore Volume', hue='type', hue_order=['mid', 'pred', 'high'])

<AxesSubplot:xlabel='main_pore', ylabel='Pore Volume'>

plt.figure(figsize=(18, 12))
plt.subplot(212)
colName = 'Pore Volume'
colNameNew = colName + r" ($\mu m^3$)"
dTmp = nb.GetMatchedOrOnMatched(allSharedPoresFeatures, common_pores, 'high', 'pred')
dTmp[colNameNew] = dTmp[colName]
sns.scatterplot(data=dTmp.drop([14]), x='pore', y=colNameNew, hue='Is Detected')
plt.title("(b) Pores detected in improved-resolution data")
plt.subplot(211)
dTmp = nb.GetMatchedOrOnMatched(allSharedPoresFeatures, common_pores, 'high', 'mid')
dTmp[colNameNew] = dTmp[colName]
sns.scatterplot(data=dTmp, x='pore', y=colNameNew, hue='Is Detected')
plt.title("(a) Pores detected in mid-resolution data")

/home/bzr0014/git/additive2/additive2/notebooks/pore_preformance.py:13: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  tmp['has_match'] = ~tmp[predCol].isna()
/home/bzr0014/git/additive2/additive2/notebooks/pore_preformance.py:13: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  tmp['has_match'] = ~tmp[predCol].isna()

Text(0.5, 1.0, '(a) Pores detected in mid-resolution data')

filt = dTmp['Is Detected']
from scipy.stats import ttest_ind

ttest_ind(dTmp.loc[filt, 'Pore Volume'], dTmp.loc[~filt, 'Pore Volume'])

Ttest_indResult(statistic=5.7529138721912005, pvalue=1.3569643216397724e-06)

from scipy.spatial.distance import pdist, cdist

a = []
b = []
c = []
for _, (i, j) in common_pores[['high', 'pred']].dropna().iterrows():
    tmp1 = pts_high[clusters_high_[i]]
    tmp2 = pts_pred[clusters_pred_[j]]
    tmp = cdist(tmp1, tmp2) == 0
    a.append(tmp.sum())
    b.append(tmp1.shape[0])
    c.append(tmp2.shape[0])
(
    np.array(a).sum() / np.array(b).sum(),
    np.array(a).sum() / np.array(c).sum(),
)   

(0.8313364055299539, 0.9504741833508957)

a = []
b = []
c = []
for _, (i, j) in common_pores[['high', 'mid']].dropna().iterrows():
    tmp1 = pts_high[clusters_high_[i]]
    tmp2 = pts_mid[clusters_mid_[j]]
    tmp = cdist(tmp1, tmp2) == 0
    a.append(tmp.sum())
    b.append(tmp1.shape[0])
    c.append(tmp2.shape[0])
(
    np.array(a).sum() / np.array(b).sum(),
    np.array(a).sum() / np.array(c).sum(),
)   

(0.6307692307692307, 0.9080068143100511)

Draft

def elongation(r, l):
    r = np.reshape(r, (-1, 1))
    l = np.reshape(l, (1, -1))
    a = 2 * np.pi * r * l
    v = np.pi * r * r * l
    return a**3/v**2

x = np.arange(1, 10)#.round(2)
res = elongation(x, x)

fig, ax = plt.subplots(figsize=(10, 10))
x_, y_ = np.meshgrid(x, x)
ax.imshow(np.log(res+1), origin='lower')
ax.set_xticks(x);
ax.set_xticklabels(x);
ax.set_yticklabels(x);
ax.set_xlabel('L')
ax.set_ylabel('R')

/home/bzr0014/miniconda3/envs/additive2/lib/python3.6/site-packages/ipykernel_launcher.py:6: UserWarning: FixedFormatter should only be used together with FixedLocator
  

Text(0, 0.5, 'R')

elongation(10, 1), elongation(1, 10), elongation(1, 1), elongation(10, 10)

(array([[2.51327412]]),
 array([[251.32741229]]),
 array([[25.13274123]]),
 array([[25.13274123]]))

plt.plot(res[:, 5])

[<matplotlib.lines.Line2D at 0x7fa3d1c2d7f0>]